; docformat = 'rst'
;
;+
;
; :Purpose:
;   Produce strong oxidant field
;
; :Inputs:
;   Sim_name: (string) Simulation name
;
; :Requires:
;   {input_directory}/EUM/grid.inp.sav
;   {input_directory}/regions.sav
;   
; :Outputs:
;   {input_directory}/EUM/OX_0001.nc
;
; :Example::
;   IDL> cels_mozart_strong_oxidant, 'SF6'
;   
; :History:
; 	Written by: Matt Rigby, MIT, Oct 12, 2011
;
;-
Pro cels_mozart_strong_oxidant, sim_name
    
  compile_opt idl2  
  
  nLevels=28 ;MLR - Assumes 28 vertical levels
  
  restore, cels_filestr(/Input, sim_name + '/regions.sav')
  restore, cels_filestr(/Input, sim_name + '/EUM/grid.inp.sav')

  LonSize=n_elements(lon)
  LatSize=n_elements(lat)

  OX=replicate(1., LonSize, LatSize, nLevels, 3)

  nRegions=n_elements(lon_local->keys())

  lon_domain_all=lon_domain_EUM
  lat_domain_all=lat_domain_EUM
  
  for ri=0, nRegions-1 do begin
    lon_domain_all[*, ri]+=[-10., 10.]  ;EXTENDED REGIONS
    lat_domain_all[*, ri]+=[-10., 10.]  ;EXTENDED REGIONS
  endfor


  for ri=0, nRegions-1 do begin
    if lon_domain_all[0, ri] gt lon_domain_all[1, ri] then begin
      whlat=where(lat ge lat_domain_all[0, ri] and lat lt lat_domain_all[1, ri])
      whlon=where(lon ge lon_domain_all[0, ri])
      OX[whlon[0]:whlon[-1], whlat[0]:whlat[-1], *, *]=0.
      whlon=where(lon le lon_domain_all[1, ri])    
      OX[whlon[0]:whlon[-1], whlat[0]:whlat[-1], *, *]=0.
    endif else begin
      whlat=where(lat ge lat_domain_all[0, ri] and lat lt lat_domain_all[1, ri])
      whlon=where(lon ge lon_domain_all[0, ri] and lon lt lon_domain_all[1, ri])
      OX[whlon[0]:whlon[-1], whlat[0]:whlat[-1], *, *]=0.
    endelse
  endfor
  
  fid=ncdf_create(cels_filestr(/Input, sim_name + '/EUM/OX_0001.nc'), /clobber)
  
  LonDID  = NCDF_DIMDEF(fID, 'lon', LonSize)
  LatDID  = NCDF_DIMDEF(fID, 'lat', LatSize)
  LevDID  = NCDF_DIMDEF(fID, 'lev', nLevels)
  TimeDID  = NCDF_DIMDEF(fID, 'time', 3)
 
  lonID = NCDF_VARDEF(fId,'lon', [LonDID], /DOUBLE)
  latID = NCDF_VARDEF(fId,'lat', [LatDID], /DOUBLE)
  levID = NCDF_VARDEF(fId,'lev', [LevDID], /DOUBLE)
  timeID = NCDF_VARDEF(fId,'time', [timeDID], /DOUBLE)
  dateID = NCDF_VARDEF(fId,'date', [timeDID], /LONG)
  datesecID = NCDF_VARDEF(fId,'datesec', [timeDID], /LONG)
  varID = NCDF_VARDEF(fId,'OX_VMR_inst', [LonDID, LatDID, LevDID, timeDID], /DOUBLE)

  NCDF_ATTPUT, fId, varID, "long_name", "Instantaneous OX VMR"
  NCDF_ATTPUT, fId, varID, "units", "VMR"

  NCDF_CONTROL, fId, /ENDEF
  
  ; Write the data
  NCDF_VARPUT, fId, lonId, lon
  NCDF_VARPUT, fId, latId, lat
  NCDF_VARPUT, fId, levId, indgen(nLevels)
  NCDF_VARPUT, fId, timeId, [julday(1, 1, 1980, 0), julday(1, 1, 2000, 0),julday(1, 1, 2020, 0)]
  NCDF_VARPUT, fId, dateID, [19800101, 20000101, 20200101]
  NCDF_VARPUT, fId, datesecID, replicate(0L, 3)

  NCDF_VARPUT, fId, varID, OX
  
  ncdf_close, fid

End
